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Abstract. The torque felt by a non-accreting protoplanet on a circular orbit embedded in a uniform surface density 
protoplanetary disk is analyzed by means of time-dependent numerical simulations. Varying the viscosity enables 
one to disentangle the Lindblad torque (which is independent of viscosity) from the corotation torque, which 
saturates at low viscosity and is unsaturated at high viscosity. The dependence of the corotation torque upon the 
viscosity and upon the width of the librating zone is compared with previous analytical expressions, and shown to 
be in agreement with those. The effect of the potential smoothing respectively on the Lindblad torque and on the 
corotation torque is investigated, and the question of whether 3D effects and their impact on the total torque sign 
and magnitude can be modeled by an adequate smoothing prescription in a 2D simulation is addressed. As a side 
result, this study shows that the total torque acting on a Neptune-sized protoplanet is positive in a sufficiently 
thin, viscous disk (H/r < 4%, a > 10 -2 ), but the inward migration time of smaller bodies is still very short, 
making it unlikely that they reach the torque reversal mass before having migrated all the way to the central 
object. 
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1. Introduction 

The formation of planets is likely to occur in the proto- 
planetary disks surrounding the very young stars. One of 
the most widely accepted formation scenario is the core 
accretion one, whereby a solid core is slowly built up by 
the accretion of solid material (rocks and ice) and subse- 
quently, if its mass becomes large enough, by the rapid 
accretion of the surrounding gas. As the core builds up, 
it gravitationally interacts with the surrounding nebula, 
and as a result of the exchange of angular momentum 
and energy with the nebula, it undergoes a change in 
its orbital parameters. The basis of the gravitational in- 
teraction between a disk and an orbiting point-like per- 
turber was worked out about twenty years ago (Goldreich 
& Tremaine 1979 and 1980, Lin & Papaloizou 1979) and 
since then refined in great detail (Ward 1986 and 1988, 
Artymowicz 1993, Ward 1997). The most striking pre- 
diction of these works was that a protoplanet embedded 
in a protoplanetary disk would undergo a rapid orbital 
decay or inward migration, and this expectation has re- 
ceived strong observational support from the discovery of 
short-period giant planets (Mayor & Qucloz 1995), the 
formation of which is unlikely to have occurred in situ. 
Most of the work done on the planet-disk tidal interac- 
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tion has been done in the linear regime. Each azimuthal 
Fourier component of the planet perturbing potential can 
be shown to contribute sizably to the disk-planet tidal in- 
teraction wherever its frequency in the local frame either 
vanishes or matches ±«, the epicyclic frequency (which co- 
incides with the Keplerian frequency in a Keplerian disk 
such as a protoplanetary disk, where self-gravity, at least 
in the central parts, in considered as unimportant). The 
angular momentum exchange rate through the terms for 
which the perturber frequency in the local frame is +n 
(resp. —k) is called the outer Lindblad torque (resp. the 
inner Lindblad torque) , while the terms corresponding to 
a vanishing perturbed frequency in the local frame corre- 
spond to the corotation torque. If one makes the assump- 
tion of a circular orbit (which is reasonable as it can be 
shown that the global effect of the resonances is an ec- 
centricity damping in most cases, see e.g. Papaloizou & 
al. 2001 and references therein), then the Outer Lindblad 
Resonances (OLR) fall exterior to the orbit, the Inner 
Lindblad Resonances (ILR) fall interior to it, and the coro- 
tation resonances nearly coincide with the orbit (they ac- 
tually lie slightly inside of this latter as the disk is slightly 
sub-Keplerian due to the partial support of gravity by a 
radial pressure gradient), and the corresponding torque is 
called for obvious reasons the co-orbital corotation torque. 
The outer Lindblad torque can be shown to be negative, 
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while the inner Lindblad torque is positive. There exists 
a mismatch between both (almost always in favor of the 
outer Lindblad torque, Ward 1997) which scales as h, the 
disk aspect ratio, and which leads to the inward migra- 
tion. The Lindblad torque has been shown to be inde- 
pendent of viscosity (Meyer- Vernet & Sicardy 1987, and 
Papaloizou & Lin 1984). This is true as long as the sur- 
face density profile itself is not sizably modified under the 
action of the one-sided Lindblad torque, i.e. as long as the 
relative depth of the dip opened around the orbit by the 
planet is small. This assumption is fulfilled in the present 
work except for the highest masses a nd the lowest vis- 
cosities (see also discussion at section 4.4). On the other 
hand, the corotation torque in the linear regime does not 
correspond to an angular momentum flux carried away 
by waves. The angular momentum rather accumulates at 
corotation, and the corresponding torque appears as a dis- 
continuity in the angular momentum flux at corotation 
(Goldreich & Tremaine 1979). Ward (1992) has shown, by 
summing the corotation resonances after an adequate ver- 
tical averaging of the perturbing potential, that the coro- 
tation torque is at most comparable to the differential 
Lindblad torque, whereas Korycansky & Pollack (1993) 
have shown by solving numerically the linearized equa- 
tions for a steady state perturbation of the disk that the 
actual solution was even smaller than the analytical es- 
timates. The corotation torque scales as the gradient of 
specific vorticity across corotation d{Ti/B)/dr, where B 
is the second Oort's constant and is 0/4 in a Keplerian 
disk. Ward (1991) has shown that the corotation torque 
comes from the exchange of angular momentum between 
the planet and the near-by fluid elements which move from 
one side to the other of their horseshoe streamline, and 
that the dependency upon the specific vorticity gradient 
was due to an uneven mapping of the fluid elements be- 
tween the outer and inner leg of their horseshoe stream- 
line. In the linear regime, i.e. as the planet mass tends to 
zero, the horseshoe zone width tends to zero, while the li- 
bration time of the fluid elements on the horseshoe stream- 
lines tends to infinity, which is in agreement with the pic- 
ture of a localized and constant discontinuity of the angu- 
lar momentum flux at corotation. In the finite mass regime 
the situation is quite different however. The libration time 
can be much shorter than the migration time across the 
horseshoe region, and also much shorter than the viscous 
diffusion time, in which case libration removes the specific 
vorticity gradient (Ward 1992), and the corotation torque 
therefore vanishes (saturates) after a few libration times. 
In a sufficiently viscous disk however, the viscous diffusion 
and the radial transport of material across the horseshoe 
region may inhibit the corotation torque saturation. A way 
of evaluating the co-orbital corotation torque as a function 
of viscosity has been given by Masset (2001), for a planet 
held on a circular orbit in a uniform surface density disk, 
and for a steady flow in the planet frame. As the material 
set in libration in a steady flow is trapped, the total torque 
on the librating fluid elements region vanishes. This torque 
is the sum of the viscous torque on the trapped region and 



the gravitational torque of the planet, which is opposite 
to one component of the corotation torque. As the vis- 
cous force is a contact action, the integral of the viscous 
torque on the trapped region reduces to the torque ex- 
erted on its boundary (the separatrix between circulating 
and librating streamlines) by the outer and inner disks. 
The second and last component of the corotation torque 
is given by the angular momentum lost by the fluid el- 
ements which participate in the global accretion of the 
disk material onto the central object (accretion of which 
are excluded only the fluid elements of the librating re- 
gion), when they flow from the outer to the inner disk 
and undergo one horseshoe like close encounter with the 
planet. The corotation torque can therefore be expressed 
(in the steady state case) only with the knowledge of the 
flow properties at the separatrix. The flow properties in 
the librating region can be split into an even and an odd 
part of the perturbed density (w.r.t the distance to coro- 
tation). The former comes from the perturbations of the 
surface density and azimuthal velocity profiles under the 
action of the one-sided Lindblad torque, while the latter 
comes from libration^} Masset (2001) gives the following 
expression for the corotation torque: 



r c = r^ + r^ 



(i) 



where F^ is the main term, which comes from the odd 
part of the perturbed density and which therefore is linked 
to libration and exhibits a dependence on the viscosity: 
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where F{z s ) is a function of the viscosity and the horse- 
shoe zone width x s (in which, contrary to Masset 2001, 
a factor 4 has been introduced for convenience, cf. in- 
fra) , while F^ is an additional term which comes from the 
even part of the perturbed density, and which introduces 
a coupling between the one-sided Lindblad torque and the 
corotation torque: 



rS = -G(x s )r LR 



(3) 



where G(x s ) = 0(1) is a function which depends on the 
exact position of the separatrix, and where IYr is the 
one-sided Lindblad torque (the arithmetical average of the 
Outer and Inner Lindblad torques absolute values). The 
fact that this additional term does not exhibit any de- 
pendence on the viscosity may seem paradoxical, as one 
expects the co-orbital corotation torque to vanish in an in- 
viscid disk. This paradox is only apparent however. This 



1 Libration will tend indeed to flatten radially across the 
horseshoe region the profile of any quantity conserved along a 
streamline. This is the case of the specific vorticity _B/E in an 
inviscid disk. Libration therefore tends to impose a decreasing 
profile of E around corotation in an inviscid disk. Only in the 
special case of the shearing sheet, in which one neglects the 
radial variation of B, does libration tend to eliminate in-out 
surface density differences. 
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term comes from the even terms in the perturbed den- 
sity (and odd in the perturbed azimuthal velocity) , which 
scale as v -1 , therefore this dependence cancels out with 
the v overall dependence of the viscous torque on the sep- 
aratrix. This results holds as long as the even perturbed 
density does not saturate (i.e. as long as the dip is shallow 
enough and does not correspond to a gap) . In the strictly 
inviscid limit (and provided that in this limit one can have 
a steady state situation), the planet (which in this artifi- 
cial situation is held on a fixed circular orbit) opens a gap, 
and therefore the even perturbed density term saturates 
and the additional term given by Eq. (g) vanishes. 

The purpose of this paper is to investigate this situa- 
tion and to check by means of numerical simulations the 
validity of Eqs. (§) and (g). For that purpose, numerical 
simulations are performed over a large range of viscosi- 
ties. Other questions addressed in this paper are the abil- 
ity of 2D numerical codes to predict correctly the total 
torque acting on an embedded point-like object through 
the choice of a correct value of the potential smoothing 
length, and the possibility of migration reversal (i.e. out- 
ward) due to the corotation torque, although linear regime 
studies seem to conclude that the (positive) corotation 
torque is too weak to counteract the (negative) differ- 
ential Lindblad torque. Although the additional term of 
Eq. (g) is too small to cancel out the differential Lindblad 
torque by itself, it is of interest to examine the situation 
for weakly embedded objects where the corotation torque 
is lift up by the additional term of Eq. (g), and where 
the differential Lindblad torque may be reduced with re- 
spect to its linear value (see e.g. Miyoshi et al. 1999, Lin 
& Papaloizou 1979). 

It should be noted that global simulations are needed 
to describe properly the corotation torque. Any local de- 
scription, such as the one provided by a box or a wedge 
centered on the planet, will prevent saturation if in- 
flow/outflow boundary conditions are used in azimuth, 
which inject "fresh" material with an unperturbed den- 
sity and velocity profile, while the use of periodic bound- 
ary conditions in these configurations artificially shortens 
the libration time. Furthermore, the use of the shearing 
sheet formalism, which usually neglects the radial deriva- 
tive of the unperturbed disk specific vorticity, leads to a 
dependence of the corotation torque on the radial deriva- 
tive of the surface density, thus switching off the corota- 
tion torque in the convenient situation where the surface 
density is uniform. 

2. Notations 

The planet orbital radius is r p , and its orbital frequency 
ftp. The position of a fluid element in the disk is repre- 
sented by its polar coordinates radius r and azimuth 8, 
counted rotation-wise, with its origin in the planet direc- 
tion. The distance of a fluid element to the planet orbit is 
x = i — r p , and its dimensionless counterpart is x = x/r p . 
The Keplerian frequency is fiif(r), and the disk mate- 
rial orbital frequency is f2(r). The disk kinematic viscos- 



ity is v and its dimensionless counterpart is v — v /(r 2 Q, p ). 
The a parameterization (Shakura & Syunyaev 1973) will 
also be used, with v = aHc s , H being the disk ver- 
tical scale length and c s being the sound speed. The 
disk surface density is denoted with E, while its unper- 
turbed uniform value is with So. The planet mass is m p , 
the central object mass is M„, and the ratio of both is 
q = m p /M*. The torque exerted on the planet is denoted 
F, and its non-dimensional counterpart is T = T/Tq where 
To = TTrpilpq 2 Y,o / h? , h being the aspect ratio. The per- 
turbed velocity has radial and azimuthal components v r 
and vg. 

3. Numerical aspects 

The code that was used is an Eulerian ZEUS-like 2D code 
based on a polar grid centered on the primary (Stone & 
Norman 1992, Nelson et al. 2000), and corotating with 
the planet (Kley 1998). The runs described here were 
performed using a modified azimuthal Courant condition 
(Masset 2000) in order to increase the timestep, and some 
of these runs were checked against the standard azimuthal 
transport procedure, and found to give almost identical 
results. The grid consists of Ng — 450 by N r = 143 
zones, uniformly spaced in azimuth and radius. The planet 
lies at radius r p — 1 and azimuth = 0. The cen- 
tral star mass and gravitational constant are respectively 
M* = 1 and G = 1, and the time unit is chosen to be 
(r p /GM»,) 3 / 2 = ilp 1 , so that the planet orbital period in 
our unit system is 2n. The grid outer boundary is chosen to 
be at r — 2.5, while the inner boundary is at r = 0.504651, 
so that the planet is located just at the center of a zone, in 
order to limit any bias in the torque and in the co-orbital 
dynamics due to an uneven placement of the planet w.r.t. 
the grid. 

3.1. Basic equations 

The equations solved by the code are: 



The continuity equation, which reads: 
<9£ 1 d{rv r Y) 1 d{v Y>) 



dt 



Or 



00 







(4) 



— The Navier-Stokes equations, which reads respectively 
for v r and vg: 



dv r dv r vg dv r v'g 

dt r dr r 89 r 

and : 

dvg dvg Vg dvg V r Vg 

dt r dr r d0 r 



1 dp dej) 
£ dr dr 



fr 



(5) 



^rd9 rd0 £ 1 ' 



where p is the vertically integrated pressure, f r and fg 
respectively the radial and azimuthal component of the 
viscous force per unit surface, and <fi is the gravitational 
potential. 
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The equation of state is that of an isothermal gas with 
sound speed c s : 



P = c s S 



(7) 



The gravitational potential includes the potential of 
the primary, the potential of the protoplanet, and an 
indirect term which arises from the fact that the non- 
rotating frame centered on the primary is not inertial: 

,. GM» Gm„ 

<t>{r,0) = 



r 

Gm 



(r 2 + r 2 
—r cos 9 



2rr p cos 9 



2)1/2 



r 



v 



GS(r', 



■ cos(<9 - 9')r'dr'd9' (8) 



disk ' 

In this expression, e is the smoothing length of the 
protoplanet potential, and otherwise stated it amounts 
to 60 % of the disk "thickness" H = c s /il. 
— The viscous force is derived from the viscous stress 
tensor as follows: 

1 dr rr f, 
r dtp r 
1 



_ 1 d(rr rr ) 

Jr 

r or 
1 9(rr 0r ) 

u = -■ 



(9) 
(10) 

or r o<p r 
where the components of the viscous stress tensor are 



r 



= 2r\D rr - g??V.v 



(11) 



— 2rjD<p<p 
where 

dv T 
dr 



3 



77 V. v 



Drr 



1 dv^ v r 

^>^>4> = - -£T H 

r 00 r 



(12) 



D r 



d_ 

dr 



1 dv r 
r da 



and r\ = T,v is the vertically integrated dynamical vis- 
cosity coefficient. 

3.2. Boundary conditions 

The boundary conditions are 2-7r-periodic in 9 in order to 
account for the disk geometry. As the viscosity can be large 
in some runs, it is important to take adequate boundary 
conditions, otherwise the fast radial redistribution of disk 
material can affect the slope of the surface density profile, 
which in turn can significantly alter the corotation torque 
magnitude. For this reason, it is important : 

— to have a source of material with surface density So a t 
large radius in order to avoid outer disk depletion; 

— to allow the outflow (inwards) of disk material at the 
inner boundary only if the disk surface density is larger 
or equal to So (otherwise a positive radial gradient of 
surface density develops, and one may overestimate the 
corotation torque in these conditions). 



3.3. Initial setup 

In all the runs presented here, the unperturbed disk sur- 
face density is uniform and is, in our unit system, So = 
10 -4 . This translates into: 



S = S ■ 8.9 • 10 6 



Vl AUJ 



g.cm 



(13) 



Note that since the results presented here are scaled by 
To, they are independent of the actual value of So (except 
for the disk potential indirect term of Eq. (|J), which for 
the value chosen above is negligible). 

The disk temperature profile is constant in time and 
is chosen such that the disk aspect ratio h = H/r be 
uniform. The initial azimuthal velocity is computed ac- 
cordingly and is slightly sub-Keplerian due to the partial 
support of gravity by the radial pressure gradient: 



GAL 



(I - h 2 



1/2 



(14) 



The initial radial velocity is set to zero. 
3.4. Initial parameters 

A number of runs have been performed varying the planet 
mass, the disk aspect ratio, the viscosity, and in some 
cases the smoothing prescription or the resolution (which 
was then chosen twice higher). Each run consisted of 
120 planet orbits (which was assumed to be sufficient to 
reach a steady state in the planet frame). For a given 
choice of the planet mass and disk aspect ratio, usually 13 
runs were performed for different values of the viscosity, 
logarithmicly spaced: 



10 -7+ 2! /7 (1 < i < 13) 



(15) 



The set of main runs is described at Tab. [l| The last 
cell of the h = 6 % line is empty: the corresponding sets of 
runs have not been performed, as the corresponding run 
would correspond to a horseshoe region radially not wider 
than one zone, and even if one disregards finite resolu- 
tion effects, the mismatch between corotation and orbit is 
not negligible compared to the horseshoe zone half-width, 
which is one of the condition of validity of Eqs. (||) and (||). 
The first cell of the h = 3 % line is also empty, as the situ- 
ation it would describe would rather correspond to type II 
migration, i.e. the planet opens a gap around its orbit and 
the corotation torque is switched off, except for the largest 
viscosities. 

Subsidiary runs have been performed in which either 
the smoothing or the resolution has been changed with 
respect to the main runs. These runs are presented at 
Tab. |. 

Subsidiary runs have been performed in which cither 
the smoothing or the resolution has been changed with re- 
spect to the main runs. The runs with a modified smooth- 
ing length are presented at Tab. || whereas the high res- 
olution runs are presented at section [| 
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Table 1. The set of 247 main runs. A run with mass 
ratio q, aspect ratio h and viscosity i>i is labeled R.m l h , 
where h is in percent and where m is the round-off of the 
planet mass, in earth masses, if M* = M@. For instance, 
run R6§ stands for the run with q — 1.67 • 10~ 5 , aspect 
ratio h — 0.03, and viscosity v%. The bulk of the runs 
are Rm l h , with m — 2,6,11,17,30, h = 3,4,5,6 % and 
i = 1 to 13. Runs with m = 30 have been performed 
for extrapolation purposes only, since for such a mass a 
protoplanet is expected to rapidly accrete the surrounding 
gas. 



<s 


3 4 5 6 


6.67- l(T b 
1.67 ■ 10~ 5 
3.33 • 10~ 5 
5 ■ 1CT 5 
10" 4 


^1^13 1^1^13 ^i-,i3 n/a 

^1^13 ^1^13 ^1->13 KL-113 
2A^13 ^1-»13 ^1 — 13 fr->13 
^1^13 ^1^13 ^1->13 KL-.13 

n/a ^i->i3 ^i-»i3 



Table 2. The 26 subsidiary runs. The smoothing length 
of the potential has been set to a smaller value (resp. 20 % 
and 40 % of the disk thickness) than its usual value (60 % 
of the disk thickness). 



Run name 


corresponding main run 


parameter changed 


S20R15^ 


R15s,i = 1 -> 13 


Smoothing: e = Q.2H 


S40R15^ 


R15s,i = 1 -> 13 


Smoothing: e = 0.4H 



4. Run results 

4.1. Generic disk response 

The disk response to the perturbing potential of a planet 
on a uniform circular orbit is illustrated at Fig. [|. The 
planet excites waves in the disk at all wavenumbers 
(Goldreich & Tremaine 1979), which corotate with the 
planet, and the superposition of which leads to a one- 
armed spiral shock (Goodman & Rafikov, 2001). This lat- 
ter is in advance w.r.t. the planet in the inner disk, which 
corresponds to the positive inner Lindblad torque, while it 
is delayed w.r.t the planet in the outer disk, which corre- 
sponds to the negative outer Lindblad torque. The balance 
between the one-sided Lindblad torque and the viscous 
torque, for a steady flow in the planet frame, enables one 
to evaluate the depth of the dip which surrounds the or- 
bit (Masset 2001), which, as long as it is shallow, scales 
proportionally to the planet mass square and inversely 
proportionally to the disk viscosity. The dynamics of the 
co-orbital region, i.e. the horseshoe region, is better seen 
in the (9, r — r p ) plane, as shown in the following sections. 



Surface density map at end of run K17 4 Surface density 




-2 -1 1 2 



Fig. 1. This plot shows the disk surface density response 
to the planet potential, at the end of run R17 4 . The rota- 
tion is anti-clockwise. The one-armed shock triggered by 
the planet in the disk is clearly visible, as well as the shal- 
low dip surrounding the orbit, the relative depth of which 
amounts to ~ 15 % of the disk surface density. 



Planet mass: 16.7 M ffl ; v; H/r: 0.05 
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Fig. 2. Streamlines aspect in the co-orbital region for 
run RI75. The thick line represents the outer downstream 
separatrix. 



4.2. An illustrative example 

Let one consider the runs R17g and R33g, which corre- 
spond respectively to: q = 5 • 10~ D , h = 0.05 and v = 
1.93 • 10~ 7 and to: q = 10~ 4 , h = 0.05 and v = 1.93 • 10" 7 . 
The aspect of the streamlines at time t = 754 SI" 1 in the 
frame corotating with the planet in either case is repre- 
sented at Fig. |^ and |^. The streamlines are found inte- 
grating the velocity field, which is bilinearly interpolated 
from the values at the interfaces between adjacent zones. 
As most of the motion in the co-orbital region can be ac- 
counted for by a linear shear v y — 2Ax where A is the 
Oort's constant, this bilinear interpolation gives satisfac- 
tory results, and in particular enables one to search for 
the precise position of the separatrix between the librat- 
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Planet mass: 33.3 M ffi ; v; H/r: 0.05 



< 1.04 




Planet masses: 16.7 M„ and 33.3 M„; H/r: 0.05 



Fig. 3. Same as Fig. || for run R335. The separatrix lies 
further from the planet orbit than in run RI75. 



ing and circulating streamlines. The separatrix position is 
defined by its distance x s to the orbit in opposition with 
the planet, i.e. at azimuth 8 = it, as indicated in Figs. || 
and |^. The value of x s is determined with a dichotomic 
search. If a guess of x s gives a circulating (resp. librat- 
ing) streamline, a smaller (resp. higher) value is tried. The 
search is iterated until a sufficient precision is achieved. 
Because the bilinear interpolation of the velocity field is 
well adapted to the kinematics of the co-orbital region, 
the precision that can be achieved on the separatrix posi- 
tion is much smaller than the radial zone size, and in that 
sense such a precision is illusory. However, this allows to 
investigate the behavior of the librating zone width when 
varying any parameter of the problem, the grid resolution 
being fixed. An estimate of the error on the corotation 
torque and on the separatrix position due to finite resolu- 
tion effect can be found in Appendix |X| For the run RI75, 
the separatrix position is found to be: x s = 0.0424 ± 10~ 4 , 
while for the run R33,!;, the separatrix position is found to 
be: x s — 0.0597 ± 10~ 4 , as can be estimated from Fig. || 
and ||. Therefore the viscous time-scale across the horse- 
shoe zone half-width is: 



r viar = -i ~ 3.1 • lO 3 ^- 1 for Rl7l 



3;/ 



6.2 ■ lO 3 ^ 1 for R33^ 



while the outermost horseshoe turnover time is: 



877 O-i 



200 ftp 1 for R17\ 
140 ftp 1 for R33 1 



(16) 
(17) 

(18) 

(19) 
(20) 



As the viscous time-scale across the horseshoe region is 
much higher than the horseshoe turnover time in either 
case, one expects the corotation torque to saturate. The 
total torque exerted by the disk on the planet as a func- 
tion of time is represented at Fig. |J. The total torque 
tends towards a constant value in time on a time-scale 




40 60 

Time (orbits) 



Fig. 4. Torque on the planet of run K17l (solid line) and 
run R30° (dashed line) as a function of time. The outer- 
most horseshoe turnover time is indicated in either case. 



The torque is normalized by Tq 



^FV/Zi 3 . The 



glitch around £ = 10 orbits corresponds to an edge ef- 
fect of the smoothing window (these results are smoothed 
over a temporal window of about 10 orbits). 




Time (orbits 



Fig. 5. Torque on the planet of runs R17| as a function of 
time, for 1 < i < 12. The steady state regime is reached 
faster as the viscosity increases, and tends towards a par- 
tially saturated state function of the viscosity. The glitch 
at t = 4 orbits is an edge effect of the smoothing win- 
dow. For this series of runs only the integration time was 
150 orbits in order to check that the total torque value 
does not vary significantly after 120 orbits. 



of the order of the outermost horseshoe turnover time. 
The initial torque value corresponds to the sum of the 
differential Lindblad and corotation torque, this later be- 
ing positive, since the surface density is uniform on the 
co-orbital region. Incidently one can notice that in both 
cases the initial corotation torque value almost cancels out 
exactly the limit value at large time, and that this limit 
value at large time is not the same for both planets. Fig|| 
shows the behavior of the total normalized torques as a 
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function of time for the case q = 5 • 10~ 5 , h = 0.05, for 
different viscosities. The bottom curves correspond to low 
viscosity, and the limit value at large time increases with 
viscosity. This limit value reaches a maximum which corre- 
sponds roughly to the value at t = 0, which corresponds to 
a totally unsaturated torque, and then decreases beyond 
the cut-off viscosity, as the topology of the co-orbital re- 
gion undergoes a change that will be discussed further (see 
Masset 2001). 

4.3. Dependence of corotation torque on the viscosity 

Fig. |^ shows the value of the normalized total torque act- 
ing on the protoplanet as a function of the reduced viscos- 
ity, evaluated from runs RI75. In each case the value of 
the torque is averaged on the time interval [90, 120] orbits, 
in order to get rid of the transitory regime of the begin- 
ning, which corresponds, as we mentioned in section 4.2, 
to the time needed for the corotation torque to (possibly 
partially) saturate. This figure also shows the best fits by 
the functional dependences given by Masset (2001), which 
are respectively: 



Planet mass: 16.7 M ffl ; H/r: 0.05 



l a — i , 



1 



c 



1 



g{z s ) 



zf zjg'{z s ) 



for the solid line, 

g{z s ) 



lC-b 



° z s g'(z s ) 
for the dashed line, and 



c 



pmax 
1 C 



g'(0) 

g'(zs) 



(21) 



(22) 



(23) 



- 1 / 3 , where TgJ ax 



for the dotted line, where z s = x s (2nis) 
(9/8)Eo^p2 ; si an d where g(z) is a linear combination of 
the Airy functions Ai and Bi which cancels out for z = 0, 
e.g.: 



g(z) = Bi(z) - v / 3Ai(z) 



(24) 



The fits presented at Fig. ^ have the following functional 
form: 



r(z) = r min + t c (z) 



(25) 



where Tc is one of the three functions given by Eqs. (pl|), 
( p2| ) or (^3|) and therefore has one free parameter, the floor 
level r m ; n , which corresponds to the case of a saturated 
corotation torque and therefore to the sum of the differen- 
tial Lindblad torque and viscosity independent additional 
term of Eq. (||) , which scales with the one-sided Lindblad 
torque, and which turns out to be r m i n = 0.050. Although 
x s is not a fit free parameter, the strong dependence of 
pmax U p Gn Xs ancl relatively high variation with az- 
imuth of the separatrix distance to the orbit enables one 
to consider it as a free parameter within a narrow range 
(e.g. here in the range 0.038 < x s < 0.043, which can be 
gathered from examination of Fig. |2|) . The value used for 
the fit of Fig. | falls within this range of values, and it 




Fig. 6. Torque on the planet of runs R17| (diamonds). 
The value of the separatrix distance used for the fits is 
x s = 0.0388. The vertical three-dot-dashed line shows the 
cut-off viscosity. Note that the fit maximum value at —0.01 
is also the torque value at t = (c/. Fig. |J), and corre- 
sponds to a fully unsaturated corotation torque. 

corresponds to the separatrix distance at 9 ~ 1.6 rad. It 
should be noted that if one discards the lowest viscosity 
points for the fit (for which residual numerical viscosity 
effects may affect the results), then one can get a good 
fit by Eq. ( pi] ) as well, with a slightly larger value of x s 
which still falls in the range 0.038 < x s < 0.043. Despite 
of the uncertainties introduced by the low viscosity points, 
the fit correctly accounts for a number of properties of the 
corotation torque main term (its absolute value, the half 
saturation for v ~ 0.05x 3 , and the fact that the partial 
saturation regime spans about two orders of magnitude of 
viscosity) . 

The behavior at high viscosities differs markedly from 
the fit by equation (^H|) . Indeed a cut-off of the corotation 
torque is expected when the viscous drift time of a fluid 
element across the outer half of the horseshoe region is 
shorter than half of its libration time. In that case the re- 
gion of fluid elements trapped in the co-orbital region has 
an azimuthal extent shorter than the azimuthal extent of 
the horseshoe region. Figs. ^ and |] represent the stream- 
lines for two cases of high viscosity (resp. v — 3.7 • 10 -5 
and v = 5.2 • 10~ 4 ). The outer downstream separatrix of 
Fig. |?] roughly coincides with that of Fig. ||. The material 
lying interior to it is not librating however, as is illus- 
trated by the dashed streamline. This material undergoes 
one horseshoe like close encounter with the planet and 
goes to the inner disk. The distance between the outer 
downstream separatrix and the boundary of the trapped 
region (given for r > r p by the outer upstream separa- 
trix) increases with the viscosity and eventually becomes 
as large as the horseshoe region. In that case, the material 
can go from the outer disk to the inner disk without any 
close encounter with the planet. Fig. |] shows the example 
of two open streamlines (dotted and dashed lines) which 
correspond to two fluid elements which reach the inner 
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Normalized torque vs. viscosity for H/r = 0.03 




Fig. 7. Streamline aspect for the run R17I, which cor- 
responds to the peak value of the total torque of Fig. || 
The four thin solid lines represent librating streamlines. 
One is almost exactly closed, the other ones not quite, 
as a strictly steady state in the planet frame has not 
been reached yet. The set of trapped fluid elements still 
have an azimuthal extent which is just smaller than 
2n. The dashed line corresponds to a unique streamline, 
which originates in the outer disk and eventually reaches 
the inner disk. The thick solid line represents the outer 
downstream separatrix between circulating and librating 
streamlines. The shaded area shows the set of librat- 
ing fluid elements, the boundary of which (the thick dot 
dashed line) is also the outer upstream separatrix and in- 
ner downstream separatrix. 



Planet mass: 16.7 M ffi ; v= 1.4x10 4 ; H/r: 0.05 




Fig. 8. Streamline aspect for the run R17g 3 . The three 
thin solid lines represent librating streamlines. They 
are approximately closed. The azimuthal extent of the 
trapped material (shaded area) is only about 2 rad, and its 
boundary roughly coincides with the outer closed stream- 
line. The dotted and dashed lines represent two open 
streamlines originating in the outer disk. The thick solid 
line represents the outer downstream separatrix. 
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Fig. 9. Total torque exerted on the planet as a function of 
viscosity, for a disk aspect ratio h = 0.03 and five different 
values of the planet mass. 



disk without any orbit crossing close encounter with the 
planet, a situation that is impossible at low viscosities. 
Note that the distance of the outer downstream separa- 
trix to the orbit, measured for 6 ~ 0.75 (instead of 9 = ir 
in order for the viscous drift not to affect sizably the mea- 
sure at these high viscosities), is in the range [0.035; 0.04] 
in Figs. ||, and ||. Therefore this value is almost constant 
over three decades of viscosity, which justifies a posteriori 
the choice of a unique value of x s for the fits of Fig. |[ 

4.4. Other runs results 

Figs. ^| to |lj show the value of the steady state normal- 
ized total torque as a function of viscosity for different 
planet masses, respectively for disk thickness h = 0.03, 
h = 0.04, h = 0.05, and h = 0.06. ft can checked on these 
figures that the viscosity at which the torque begins to de- 
crease increases with the planet mass. This is in agreement 
with the fact that the cut-off at high viscosity occurs for 
v ~ (l/47r)xg, and that x s increases with the planet mass. 
It can also be checked that the torque saturation occurs 
at lower viscosities for smaller masses, as the saturation 
occurs for v ~ S? s . 

These results can be used to try identify the additional 
corotation torque term of Eq. (||). The main problem is to 
evaluate the residual torque, i.e. to correct the measured 
torque from the main corotation torque term. This lat- 
ter involves the separatrix distance to the orbit x s to the 
fourth power, and it is therefore very sensitive to any error 
on the estimate of this separatrix distance. Two solutions 
can be adopted: the first one consists in considering that 
the measured torque value at low viscosity is the sum of 
the differential Lindblad torque and the additional corota- 
tion torque term (in which case there is no need to evaluate 
the main corotation torque term) and the second one con- 
sists in correcting the torque measured at large viscosity 
-although much before the cut-off- from this main term 
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Normalized torque vs. viscosity for H/r = 0.04 Planet mass: 6.67 M s ; i/ = 5.2*10~ ; H/r: 0.03 




Reduced viscosity v 




Fig. 10. Total torque exerted on the planet as a function 
of viscosity, for a disk aspect ratio h — 0.04 and five dif- 
ferent values of the planet mass. 

Normalized torque vs. viscosity for H/r = 0.05 




Reduced viscosity y 

Fig. 11. Total torque exerted on the planet as a function 
of viscosity, for a disk aspect ratio h ~ 0.05 and five dif- 
ferent values of the planet mass. 

Normalized torque vs. viscosity for H/r = 0.06 
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Fig. 12. Total torque exerted on the planet as a function 
of viscosity, for a disk aspect ratio h — 0.06 and four 
different values of the planet mass. 



Fig. 13. Separatrix shape for a low mass, low aspect ra- 
tio case. The distance of the separatrix to the orbit varies 
by about 23 % over the azimuth and therefore its fourth 
power varies by a factor 2.3 with azimuth, xf and x~ 
stand for the distance to the orbit of the outer and in- 
ner separatrices, respectively. They are functions of the 
azimuth 9. In this example A9 — 1.5 rad. 

and attributing the residual to the sum of the differen- 
tial Lindblad term and the additional corotation torque 
term. This latter solution turns out to be easier and more 
reliable, despite of the difficulty in measuring x a . Indeed 
the disk response at low viscosity may involve a strongly 
non-linear disk response (the dip around the orbit tends 
to become a fully qualified gap) , which modifies the differ- 
ential Lindblad torque in a way that is not trivial to cor- 
rect. Furthermore, finite numerical viscosity effects (which 
are difficult to quantify) and the long viscous time-scale 
are two additional factors which may lead to a significant 
difference between the ideal low viscosity steady state sit- 
uation and the run results after 120 orbits. 

Evaluating the corotation torque main term involves 
a precise measurement of the value of x s . It can be seen 
in Fig. |l3] how the distance between the separatrix and 
the orbit varies with azimuth. The value for x s is chosen 
to be the arithmetical average of the outer downstream 
separatrix position at a given angle A9 from outer con- 
junction and the inner downstream separatrix position at 
the corresponding angle 27r — A9 from inner conjunction: 



= hx+(M)+xj(2n-M)] 



(26) 



The value of AO is then varied between 1 and tt rad, 
which leads to a dispersion in the corotation torque main 
term evaluation, and this leads therefore to a dispersion 
of the residual torque estimate. The corresponding resid- 
ual torque estimates are shown for the four different disk 
thicknesses in Figs. |l^ to [l7|. 

The results of the linear regression fits performed on 
these figures are shown in Tab. |^. A number of comments 
can be made from these figures and from the table: 



10 



S -0.030 - 



F.S. Masset: The co-orbital corotation torque in a viscous disk: numerical simulations 

H/r = 3 % H/r = 5 % 




Fig. 14. Residual normalized torque as a function of the 
separatrix distance, for the 3 % aspect ratio disk. The dot- 
ted line shows the linear regression fit. The error bars are 
obtained by taking respectively the maximum and mini- 
mum value of the average of xf(A9) and xj(2ir — A8), 
when A varies. 



H/r 




0.030 



Fig. 16. Residual normalized torque as a function of the 
separatrix distance, for the 5 % aspect ratio disk. The dot- 
ted line shows the linear regression fit. The white square 
shows a measurement for which the separatrix distance 
to the orbit is smaller than the radial zone width, which 
means that the librating region is narrower than two zones 
in radius. It is therefore discarded for the linear regression 
fit. 



H/r = 6 




-0.0501 

0.00 



Fig. 15. Residual normalized torque as a function of the 
separatrix distance, for the 4 % aspect ratio disk. The 
dotted line shows the linear regression fit. 



The correction is performed at the peak value of the 
total term, while the separatrix position is searched for 
at a slightly lower viscosity. The correction takes into 
account a possibly partial saturation of the corotation 
torque at its peak value. 

The function which is used to correct the main corota- 
tion term from partial saturation is the function given 
by Eq. . Inspection of Fig. || shows that using either 
Eq. ( 22 ) or Eq. ( |2^ ) would not change significantly the 
results, given the uncertainty on x s , and given the fact 
that at its peak value the corotation torque is almost 
fully unsaturated. 

The value at x s — of the linear regression fit corre- 
sponds roughly to the limit value that can be deduced 
by eye from the limit value at v — in Fig. and [HI 




Fig. 17. Residual normalized torque as a function of the 
separatrix distance, for the 6 % aspect ratio disk. The 
dotted line shows the linear regression fit. 



In the case of h = 0.05, the uncertainties are consid- 
erable, and in the h — 0.06 case, no curve reaches the 
fully saturated regime at low v. 

The mean slope Q distribution is compatible with the 
expectation. It should be mentioned however that the 
function Q{x s ) is not expected to be a constant, and 
that its value, which is such that: 1 < Q(x s ) < 7/3, 
depends on the location x s at which the separatrix 
samples the axisymmetric profiles perturbations. 
The variation with mass of the peak value of the to- 
tal normalized torque for a given aspect ratio can be 
accounted for by the slope of the fit, i.e. this variation 
is roughly equal to t7(i!. nax — £" nn )rLR. This therefore 



F.S. Masset: The co-orbital corotation torque in a viscous disk: numerical simulations 



accounts for the increasing peak value with the planet 
mass, whereas one would expect a constant maximum 
in the linear regime (x A s cx q 2 ) and a decreasing peak 
value beyond (x* oc q 4 ^ 3 ). 
6. Not all the behavior observed is entirely due to the 
additional corotation torque term. Indeed there is 
no reason to assume that the normalized differential 
Lindblad torque is a constant for a given aspect ratio. 
Although this is certainly true in the linear regime, 
it is no longer true for higher masses. Miyoshi et al. 
1999 note that in an infinitesimally thin disk, in the 
non-linear regime (the threshold of which they find at 
Rh = H/2, where Rh is the planet Hill radius), the 
one-sided Lindblad torque is smaller than its linear es- 
timate. They mention that part of the cut-off comes 
from material inside the horseshoe region, which ex- 
erts an opposite torque on the perturber, and weak- 
ens the linear estimate of the Lindblad torque. As the 
material in the horseshoe region is librating in the 
planet frame, this material in average corotates with 
the planet, and therefore participates in the corota- 
tion torque and no in the Lindblad torque. If one de- 
fines the Lindblad torque as the torque exerted by 
the circulating material on the planet, which is, in 
the steady state regime, the only definition which en- 
sures that the total torque exerted on the planet is 
the sum of the corotation and Lindblad torques (as 
the co-orbital dynamics partitions the disk into librat- 
ing and circulating fluid elements) then it is not clear 
whether the behavior observed by Miyoshi et al. is due 
to a Lindblad torque non-linear cut-off or to an avatar 
of the additional corotation torque term. If one per- 
forms however a local scattering calculation (Lin & 
Papaloizou 1979, Papaloizou & Lin 1984) over the cir- 
culating fluid elements, then one expects a q^ 1 cut-off 
of the one-sided Lindblad torque when the separatrix 
lies further than (2/3)H from the orbit, i.e. when the 
horseshoe region invades what would be otherwise the 
place of disk-perturber angular momentum exchange 
through Lindblad torques. It should be noted that on 
the present data set, the conditions x s > (2/3)H and 
Rh > (l/2)-ff are almost strictly equivalent. In the 
h = 3% case (Fig. |lj), all points but the bottom one 
correspond to situation beyond the non-linear thresh- 
old, while in the h = 6% case (Fig. [l7]), all points cor- 
respond to a situation where the one-sided Lindblad 
torque is still in the linear regime (the Hill radius of 
the most massive planet is precisely 0.03r p ). Therefore 
the behavior observed for the thinner disks is likely 
to be contaminated by a differential Lindblad torque 
non-linear cut-off (which also conspires to lift the nor- 
malized curves as the planet mass increases), whereas 
the behavior observed for the h — 0.06 disk is exempt 
of this contamination and should be due entirely to the 
additional corotation torque term of Eq. (^). 
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Table 3. Linear fit regression values for Figs. |lJto|l7]. The 

slope Q is defined by T(x s ) — T^^o + £/i s rLR, where the 
normalized value of T^R for all the disks is 0.21. 



h{%) 


slope Q 


r* s ->o 


correlation 


3 


1.61 


-0.026 


0.978 


4 


1.61 


-0.046 


0.90 


5 


1.89 


-0.060 


0.44 


6 


2.15 


-0.077 


0.997 



which scales roughly as x s Flr, in agreement with the ex- 
pectation that G(x s ) = 0(1), and Q(x s ) < 7/3. 

5. Smoothing issues and 3D expectations 

The above simulations are performed in 2D, and the disk 
finite thickness is taken into account through the use of a 
smoothing coefficient e in the potential, where e is a sizable 
fraction of the disk vertical scale length H = r—. One can 
wonder how accurate this description is, and what fraction 
T] = e/H of the disk thickness should be used in order to 
get as close as possible a result to the 3D expectations. 
Moreover, there is no reason that an adequate value of e 
for the Lindblad torque (in the sense that this value of 
e will provide a value for the one-sided and/or differen- 
tial Lindblad torque in agreement with results obtained 
with three dimensional calculations, see e.g. Miyoshi et 
al. 1999) will also be an adequate value for the corotation 
torque. The following discussion is therefore divided in two 
parts: the effect of smoothing on the Lindblad torque is 
first investigated, then secondly on the corotation torque, 
and lastly it is discussed whether Lindblad and corota- 
tion torques can be described appropriately by the same 
smoothing coefficient. 

5.1. Effect of smoothing on the Lindblad torques 

Fig. [l8] shows the ratio of the one-sided Lindblad torque 
for a smoothed potential and the one-sided Lindblad 
torque for an unsmoothed potential, as a function of the 
ratio of the smoothing length to the disk thickness. The 
one-sided Lindblad torque is evaluated by the arithmeti- 
cal average of the outer and inner Lindblad torque expres- 
sions, which are: 

r OLR = ]T rp (27) 

m=0 

and: 

+oo 

riLR = ]T r m (28) 

m=0 

where: 



Despite of the difficulty of these measures, these re- 
sults confirm the existence of a term in the total torque 



(29) 
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Fig. 18. One-sided Lindblad torque cut-off as a function 
of the relative smoothing e/H . The long dashed line rep- 
resents the approximate asymptotic solution for h — > 
deduced from Eq. (Ba). 



in which e — ±1 (for Outer/Inner Lindblad resonance 
respectively), and where: 



{rf r G e m {r)) re +2m 



2n' 71 -QpvT 3 7? 



^ ml m / 



(l + ^/i 2 ) 1 ^ 



-(30) 



in which G'„(r) is a generalized Laplace coefficient (i.e. a 
Laplace coefficient for which a smoothing is introduced): 



G£,(r) 



cos(m0) 



-de 



ir J (1 - 2rcos6> + r 2 + e 2 ) 1 / 2 ' 
and where: 

n^, = fi P [(l - h 2 ) 1 ' 2 + em-\l + m 2 ^) 1 ' 2 ]- 1 



(32) 



is the Keplerian frequency at the effective location of the 
Lindblad resonances, and: 



yip 

and we use: 



-2/3 



(r^j = ~3meil £ m yi p {l + m 2 h 2 fl 2 



(33) 



(34) 



This calculation is identical to the calculation in Ward 
(1997) except for the introduction of the smoothing of the 
perturbing potential. 

The cut-off function is shown for three different aspect 
ratios, and found to be only weakly dependent on h. At 
low aspect ratio, the functional dependence of the torque 
cut-off tends towards a fixed function, the graph of which 
should roughly coincide with the solid curve of Fig. |f. An 
approximation of the limit cut-off function when h — ► 
can be obtained by developing Eq. (|3l]) to first order in 
mh, and approximating the G e m coefficients with a stan- 
dard technique as a function of the Bessel K and Ki 
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Fig. 19. Differential Lindblad torque cut-off as a function 
of the relative smoothing e/H . 



functions. The summation over m can then be approxi- 
mated as an integral. Using the variables £ = mh and 
r\ = e/H, one can finally write the following approximate 
formula for the one-sided Lindblad torque as a function of 
the relative smoothing 7/ only: 



r (1 



(1+ £2)1/2 



4£ 2 



F 



2„2 



(i + O + O? 



(31) F i x ) 



where if is a constant and: 
1 



Ki(V5)+Ko(>/5) 



d£(35) 



(36) 



where Ko and Ki are the modified Bessel functions. The 
long dashed line of Fig. |l8| is the graph of the cut-off 
7lr(?7)/7lr(0). Fig. |l^ shows the behavior of the differ- 
ential Lindblad torque cut-off for the same set of disk 
thicknesses. This cut-off exhibits quite similar a behav- 
ior to the one-sided Lindblad torque cut-off. On Figs, [lj 
and [l9| the dot-dashed line shows the value of the smooth- 
ing for which the torque cut-off amounts to 0.43, which 
is the value quoted by Miyoshi et al. 1999 for the ratio 
of the Lindblad torque in a vertically resolved disk and 
in an infinitesimally thin disk, in the linear regime. The 
corresponding relative smoothing is 7/ = 0.76, both on the 
one-sided Lindblad torque and on the differential Lindblad 
torque. Therefore a potential smoothing length of 76 % of 
the disk thickness should be used in 2D simulations in or- 
der to give correct results for the (differential) Lindblad 
torque. The fact that this value is independent of the disk 
aspect ratio and planet mass (at least in the linear regime) 
comes from the fact that most of the torque between the 
disk and the planet comes from a zone which is at a dis- 
tance ±H from the orbit (the Lindblad resonances pile 
up for m — > oo at a distance ±§-ff from the orbit). Such 
an argument cannot be used for the co-orbital corotation 
torque which needs a different treatment. 
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Planet mass: 16.7 M„ 
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Fig. 20. Total torque acting on the planet, as a function Fig. 21. Corotation torque main term estimate as a func- 
of the reduced viscosity, for runs S20R17|, S40R17| and tion of the separatrix position. The dotted line correspond 



to a fourth power of x s , and the dashed one to a third 
power of ir s . 



5.2. Effect of smoothing on the corotation torque 

The corotation torque comes from the exchange of angular 
momentum between the planet and nearby, orbit crossing 
fluid elements (which can be either fluid elements librat- 
ing on closed, horseshoe like streamlines, or fluid elements 
participating in the global viscous accretion of the disk, 
and which pass by the planet as they go from the outer 
disk to the inner disk). As a fluid element participating in 
the corotation torque, with impact parameter x = r — r pi 
will be located at a distance — x from the orbit after a 
back-scattering by the planet, the specific angular momen- 
tum that it gives to the planet is ABxr p , and therefore this 
value does not depend on smoothing. The total value of the 
corotation torque does depend on it however, since the ra- 
dial width of the libration region depends on the smooth- 
ing: the smaller the smoothing, the further the separatri- 
ces lie from the orbit, and the larger the corotation torque. 
This behavior can be checked at Fig. [2(]. The runs for a 
planet of mass m p = 16.7 have been performed with 
three different smoothing values. At low viscosity (i.e. for 
a saturated corotation torque) the total torque is negative 
and has a larger absolute value for a smaller smoothing. 
This corresponds to expectations since the torque in this 
region (whether it includes the Lindblad torque coupling 
term of the corotation torque or not) scales as the (one- 
sided or differential) Lindblad torque. On the other hand, 
the difference between the maximum value of the torque 
with the minimum value at low viscosity should roughly 
correspond to the maximum value of the corotation torque 
main term, which scales as xf. At Fig. [2l] we plot this dif- 
ference as a function of the separatrix distance to the orbit 
x s . Despite the small number of measurements, the results 
are in good agreement with a x^ dependency of the coro- 
tation torque. This also confirms the argument exposed 
above that the corotation torque depends on the smooth- 
ing only through the value of x s . 



It is worth emphasizing the dramatic dependence of 
the corotation torque on smoothing. In Fig. ^0|, for a 
smoothing length that is 60 % of the disk thickness, the 
torque is always negative, although it can reach a value one 
order of magnitude smaller than the differential Lindblad 
torque linear estimate for an unsmoothed potential. On 
the other hand, a smoothing length that is 20 % of the 
disk thickness leads to a torque reversal over a significant 
range of viscosity (1.5 • 10~ 5 < v < 6 • 10 -5 , which corre- 
sponds to 6 • 10" 3 < a < 0.024). 

An idea of the co-orbital dynamics in a three dimen- 
sional situation can be obtained simply if one assumes 
that the motion is purely horizontal. In that case the mo- 
tion in the slice of altitude z and infinitely small thickness 
dz is equivalent to a 2D situation, in which the potential 
smoothing length is z. Furthermore, as for \z/r p \ <C 1 
the specific angular momentum radial gradient is still 
2Br = Q.Kr/2, the torque contribution of the slice is the 
same as the torque exerted by an infinitcsimally thin disk 
(with a potential smoothed over a length z) of surface den- 
sity p(z)dz. Let X s (e,(7, h) be the separatrix distance to 
the orbit in an infinitesimally thin disk in the equatorial 
plane, in which the potential is smoothed over a length e, 
the sound speed is c s = vj^h, and the planet to primary 
mass ratio is q. The fully unsaturated corotation torque 
main term in a thick disk can therefore be written as: 



9 



^lp{z)Xt 



z,q, ■ 



c s {z) 
v K 



dz 



(37) 



As this analysis is restricted to the case of isothermal 
Keplerian disks, Eq. (B7h can be recast as: 



\nlp e-^/ H ^Xj(z, q ,h)dz 



(38) 
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A H/r = 0.06 
n H/r = 0.04 
o H/r = 0.02 



0.6 o.e 

R H H 



Fig. 22. Separatrix distance to the orbit as a function 
of the relative smoothing r\ — e/H, for different as- 
pect ratios and planet masses (q = 1.67 • 10~ 5 : dashed 
lines, q = 3.33 • 10~ 5 : dotted lines, and q = 5 • 1CP 5 : solid 
lines). As expected, x s decreases as the smoothing in- 
creases, decreases as the aspect ratio increases (everything 
else being fixed), and increases as the mass increases (ev- 
erything else being fixed). 



where op 

dl)< 



Eq 



is the disk density in the equatorial plane, 
can also be written as: 



9. 



in which x^ is defined as: 
"+°° p(z)X*(z,q,h) 



So 



dz 



which can be recast as: 



e -i(z/H) 2 X f(z,q, h)dz 



(39) 



(40) 



(41) 



Therefore, if one wants a 2D simulation to give the cor- 
rect amplitude for the fully unsaturated corotation torque 
main term, one needs to use the smoothing length which 
endows the co-orbital motion with a separatrix position 
x s , where x s is given by Eq. (ffl|). The integrand of the 
right hand side contains the function X s , which has to be 
tabulated a priori with a set of runs with different smooth- 
ing lengths, for a given planet mass and disk aspect ratio. 
The same tabulation is used to infer the correct value for 
e once x s is known. Fig. shows such tabulations for 
different planet masses and disk aspect ratios. They are 
measured with a dichotomic search of the separatrix posi- 
tion on the results of low viscosity runs {p = 10~ 6 ), over 
150 orbits. Fig. p3| shows the results of this integration for 
the nine situations of Fig. [22[ It can be seen that on the 
contrary to the case of the Lindblad torque, the value of 
the correct smoothing depends on Rh/H even when this 
value is below 1/2. Furthermore, in this case, which cor- 
responds to the linear regime for the Lindblad torque, the 



Fig. 23. Corotation torque adapted smoothing as a func- 
tion of Rh/H, inferred from the tabulations of Fig. |||. 

correct smoothing length is clearly smaller (50-60 % of the 
disk thickness) than the correct smoothing length for the 
Lindblad torque (76 % of the disk thickness). The choice of 
any value between 60 and 76% of the disk thickness will 
therefore underestimate the corotation torque and over- 
estimate the Lindblad torque, and therefore in any case 
will underestimate the total torque, since the differential 
Lindblad torque is negative and the corotation torque is 
positive. Furthermore, as the corotation torque is a very 
sensitive function of the smoothing, there is little hope of 
finding a correct smoothing prescription which predicts a 
correct total torque value in a 2D simulation, especially 
in the cases of interest where the corotation torque and 
differential Lindblad torque almost cancel out. One can 
however make conservative assumptions to infer the di- 
rection of migration of a protoplanet by choosing a lower 
or upper limit for the smoothing coefficient. This is the 
object of the next section which presents four series of 
selected high resolution runs. 

6. Higher resolution runs 

As the results of the h = 0.03 and h = 0.04 runs suggested 
that it might be possible to reverse migration over a cer- 
tain viscosity range for the most massive non-accreting 
planets (10 M® and 17 M®), the corresponding runs 
have been repeated with slightly different values of the 
smoothing following the discussion of section and with 
a higher resolution in order to rule out possible finite res- 
olution effects, which could particularly affect the differ- 
ential Lindblad torque in the thinnest disks. Namely, the 
setups detailed at Tab. || have been run. 

Each setup has been run for the five values of the vis- 
cosity indicated at Tab. [|, ranging from v — 5 ■ 10~ 6 to 
v = 10~ 4 , over 60 orbits only. This amount of time may 
seem a bit short, but a look at Fig. || shows that for these 
relatively high values of the viscosity, the limit torque 
value at large time is reached on a shorter time-scale than 
the outermost horseshoe turnover time. This turns out to 
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Table 4. The series of runs listed below all have N r — 286 
and Ng — 900, the same i? max as the previous runs while 
-Rmin has been slightly modified and set to 0.509280742 in 
order for the planet to be at the center of a zone. The set 



of reduced viscosities is in each case: v = 5 • 10 



-6 



20i /4 



i £ [0,4]. The last column shows the adopted ratio r\ of 
the smoothing length to the disk thickness. 



Run name 


Aspect ratio 


Planet mass 


Rh/H 


V 


HRllg 


0.03 


3.33 • 10" 5 


0.74 


0.76 


HR11| 


0.04 


3.33 ■ 10~ 5 


0.56 


0.55 


HR17 3 


0.03 


5 • 10~ b 


0.85 


0.76 


HR17| 


0.04 


5 ■ 10~ b 


0.64 


0.55 



-0.010 f- 




H/r = 0.03 A q = 5xl0~ 5 (17 U„' 
H/r = 0.04 <> q = 3. 33xl0" 6 (10 M„) 

Reduced viscosity v 



Fig. 24. 

Normalized total torque acting as a function of viscosity for 
runs HRII3, HRII4, HR17| and HR17|. 



be also the case with these high resolution runs. The value 
for the smoothing length has been chosen according to the 
value of the ratio of the Hill sphere radius to the disk thick- 
ness Rr/H and to the results of Fig. For the largest 
ratios, a smoothing length of 76% of the disk thickness has 
been chosen. This value is close to the correct one as far 
as the corotation torque is concerned, and it is assumed 
that it is also the correct one for the Lindblad torque, 
although such values of Rh/H do not correspond to the 
linear regime for which the canonical smoothing value was 
obtained. One the other hand, for the two smallest Rh/H 
ratios, a value r\ = 0.55 was adopted. This value should 
lead to correct results as far as the corotation torque is 
concerned, and should lead to an overestimation of the 
(negative) differential Lindblad torque, therefore they are 
conservative in the sense that they give a lower limit on 
the total torque peak value. Fig. ^4] shows the normalized 
total torque at large time as a function of viscosity for 
the four series of runs. For the runs HR17 3 , HRII3 and 
HR17I, the total torque reaches a positive peak value, cor- 
responding to an outwards migration. The limit a values 
of the interval over which the migration is outwards are 



respectively: 0.01 < a < 0.1, 7 • 10~ 3 < a < 4 • 10~ 2 , 
and 8 • 10~ 3 < a < 0.037. In this last case it should be 
remembered that the choice of the smoothing coefficient 
was a conservative one, and therefore it is likely that the 
actual viscosity interval over which migration reverses is 
larger. The case of the run series HR11| shows that at its 
peak value, around a ~ 0.01, the torque value is 22 times 
smaller than the linear differential Lindblad torque esti- 
mate for an unsmoothed potential (and this again is a 
conservative estimate). 

It is noteworthy that the results of runs HR17 3 and 
HRII3, despite of the small aspect ratio, do not differ 
much from the results obtained at lower resolution. This is 
an indication that the corotation and differential Lindblad 
torques were already satisfactorily described at that reso- 
lution. 

These results indicate that the migration of Neptune- 
sized planets (as well as slightly smaller objects) is re- 
versed in thin (h < 0.04), viscous disks (with viscosities 
typically between 0.01 < a < 0.04, although it should be 
remembered that the exact extent of this range depends 
on the planet mass and disk aspect ratio). 

7. Discussion 

The migration regime corresponding to the range of 
masses, aspect ratios and viscosities investigated in this 
work corresponds mostly to type I migration regime, ex- 
cept the low aspect ratio, low viscosity and large mass 
runs, in which quite a significant gap is opened in the co- 
orbital region. As type I migration is known to be too 
fast (i.e. the migration time of a protoplanetary core is 
shorter than its build-up time), it is of interest to inves- 
tigate the migration time from these results, as a func- 
tion of aspect ratio, planet mass and viscosity. Figs. |25| 
and ^6] show the migration time, in years, of a protocore 
at 1 AU embedded respectively in a 3% and 4% aspect ra- 
tio disk. In either case the disk surface density is chosen to 
be that of the minimum mass solar nebula (1700 g.cm~ 2 , 
i.e. So = 1.9-10 , see Hayashi et al. 1985). As the results 
of the high resolution runs of section || dealt only with re- 
stricted mass and viscosity ranges, the results shown at 
Figs. ^B] and ^6] correspond to the low resolution runs (i.e. 
to Figs. H and |l^), for which the parameter coverage is 
much larger. As the results for the h = 0.04 case in the 
high and low resolution case differ sensibly, the contours of 
Fig. (and of Fig. ^Bj) should not be taken too literally. 
They do however have the merit to show the behavior 
of the migration time as a function of mass and viscos- 
ity (even if the boundary of the domain of torque reversal 
should not be taken literally) and they also have the merit 
to give a correct order of magnitude of this migration time, 
defined as: 



2r 



(42) 



This migration time has been evaluated at r p = 1 AU, 
where the protoplanetary disk is likely to be subject to the 
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Fig. 25. Migration time in years at 1 AU for a h — 0.03 
disk, as a function of planet mass and disk viscosity. The 
dotted contours, labeled with negative values, correspond 
to a positive torque and therefore to an outwards migra- 
tion, in which case the label indicates the semi-major axis 
doubling time. The dashed line is the set of positions where 
the unperturbed viscous drift rate — 3v/2r has the same 
absolute value as the inferred migration rate. Therefore 
the migration rate estimates are reliable only to the right 
of the dashed line, i.e. wherever the material flow rate 
across the orbit is mostly accounted for by the viscous 
drift. The solid thick line corresponds to a vanishing total 
torque (limit of migration reversal). The contours close to 
this limit correspond to migration times which tend to in- 
finity, and they have been blanked for numerical precision 
reasons. 
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Fig. 26. Migration time in years at 1 AU for a h = 0.04 
disk, as a function of planet mass and disk viscosity. The 
contour line style obeys the same conventions as in Fig. E3. 



1999 and references therein), it is still much shorter than 
the disk lifetime and the core build-up time. It is therefore 
very unlikely that a unique protocore reaches the torque 
reversal mass before having migrated all the way to the 
central object. One reason is that the torque reversal oc- 
curs in very thin disks, for which the absolute value of 
the differential Lindblad torque is large (it scales as ft -2 ), 
and the associated migration time-scale is short. One can 
see that the most favorable situation in both cases is for 
a = 10~ 2 (i.e. for the corresponding viscosity a given mass 
protocore has the largest migration time). For this value 
of a one gets bottleneck values of the migration time re- 
spectively 7— 8T0 4 years (for h = 3 %) and 6— 7T0 4 years 
(for ft = 4%). 

It should also be noted that, as the corotation torque 
scales with the slope of the specific vorticity gradient, the 
torque reversal that was found here for very thin disks 
is not likely to subsist if one takes a significant negative 
surface density gradient (i.e. if £ oc r~ q with q ~ 1). 

On the other hand if the specific vorticity gradient is 
even steeper than the one envisaged here, then the corota- 
tion torque, depending on how steep this gradient is, may 
well unconditionally dominate the differential Lindblad 
torque. This might be the case in the very central parts 
of the disk, where the outer disk solution has to be con- 
nected to an inner cavity (e.g. magnetically cleared). If 
the transition region is larger than the co-orbital zone of 
an infalling body, then the corotation torque acting on 
it might be able to counteract the differential Lindblad 
torque, providing another way of stopping the inward mi- 
grating bodies, as the disk there is very likely to be viscous 
and therefore very likely to be able to sustain an unsatu- 
rated corotation torque. 

As the analysis presented here is restricted to a planet 
held on a fixed circular orbit, it should be realized that 
the migration time-scale estimated above is only valid 
whenever the viscous drift rate of material across the 
orbit is much larger than the planet inferred drift rate 
2T/(M p r p Ctp), otherwise the corotation torque expression 
has to explicitly include a term proportional to d, the mi- 
gration rate, with a delay which scales as the outermost 
horseshoe libration time. The domain for which the migra- 
tion rate is negligible compared to the viscous drift radial 
velocity is shown in Figs. ^ and |2^. Therefore one can see 
that at very low viscosities [a < (1 — 3) • 10~ 3 ] the expres- 
sion for the corotation torque needs to be reconsidered, 
as most of the specific vorticity drift across the co-orbital 
region is not accounted for by the viscous accretion but 
by the migration itself. In these circumstances the total 
torque felt by the planet depends on its migration rate. 
The analyze and consequences of this feed-back will be 
presented elsewhere. 



magneto-rotational instability (Balbus & Hawley 1991), 
which endows it with a source of significant viscosity, 
which could be in the range of the values of a for which 
migration can reverse. Although this migration time is a 
few times larger than previous results (e.g. Miyoshi et al. 



8. Summary 

This paper examines the torque exerted on a protoplanet 
held on a fixed circular orbit by a uniform surface density 
and uniform aspect ratio viscous protoplanetary disk, by 
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Outer separatrix 




means of numerical simulations. A number of runs have 
been performed, varying the disk aspect ratio, the planet 

mass and the disk viscosity, which enables one to disen- i M 
tangle the differential Lindblad torque/corotation torque 

additional term from the co-orbital corotation torque main t=1 
term. The differential torque functional dependence upon 
the width of the librating fluid elements region and upon 0rbii 
viscosity is in good agreement with previous expectations. 
The behavior of the total torque as a function of viscosity 
and planet mass is presented. It is found that the curves 
of the normalized total torque as a function of viscos- 
ity are lift up as the planet mass increases, eventually 
leading to a positive torque peak value in the thinnest 
disks. This behavior is explained as due to the non- linear 
Lindblad torque cut-off and to the corotation torque addi- 
tional term, which both conspire in lifting the normalized Fig. A.l. Sketch of a shearing sheet case. A radial set 
torque as the planet mass increases. Smoothing issues are f zones is represented. The separatrices are assumed to 
discussed. The smoothing length is shown to have a differ- be horizontal and are symmetric w.r.t the orbit. The az- 
ent impact on the Lindblad torque and on the corotation imuthal velocity is set on a staggered mesh (centered in 
torque, and the corotation torque is found to depend dra- radius, staggered in azimuth), and is assumed in this sim- 
matically on the smoothing length value, which shows that plified example to be independent of azimuth, 
there is no such thing as a "magic value" for the smoothing 

length which would give unconditionally correct results m, , , , r , 

, ,. .iii. \ -L ne c °de used here enforces angular momentum conser- 

l.e. compatible with fully three dimensional calculations . n , jn , , r 

\ .... , , . r , , . vation (to the computer accuracy) hence the balance of 

A method is given to evaluate the value of the smoothing , n , , in , . , a . , 

, r , • angular momentum gamed/lost by the orbit crossing timd 

length which gives correct results for the corotation torque n , , , , , , 

. „_ . . n T f elements on horseshoe streamlines can be used as an es- 
(l.e. identical to what one expects m a 3D situation). If ,. r ,, .. m , , 

\ . . . , timatc ot the corotation torque. Ihe angular momentum 

the motion m the disk m the co-orbital region is purely . , . n , T n . , .. . . „ 

, , . ,,11,.. • • i is a zone centered variable. In this simplified situation the 

horizontal, this method should give m principle an exact , ., „ , , . , , , , , , „, 

, 1VT ' , . . ° . , . , velocity held is assumed to be the unperturbed one. Ihe 

result. Note however that the functional dependence upon ,. , . , . , , , . , . , , ,., 

. . zone radial index is denoted i, and vanishes at the orbit 

viscosity of the 3D corotation torque and the 2D one with , , , . , , . r 

J . , (as stated m section 3L the planet lies m the center of a 

the appropriate smoothing length may be different, and , ™ . , r rr r n i i i i 

, i,i ii. i r , r zone). Ihe index of the outermost zone fully embedded 

that the method exposed here is exact only for the case of . , , , ., , ... , . T , . , , , 

... in the co-orbital region is denoted im- It is assumed that 

a fully unsaturated torque, i.e. at large viscosity although , .. , .... , ,, , 

, J rr\ tt- i i • i • tne co-orbital region is wide enough so that such a zone 

before the cut-on . High resolution runs arc presented m • , m, • r ■ • , i r 

... , i-ii.i -i i exists. Ihe expression lor Im is therefore: 

which the smoothing length is chosen conservatively, and 

these runs show that Neptune-sized planets undergo a pos- . _ „ / x s 1 

itive torque in thin, viscous disks with a shallow surface %M \ A 2 

density profile, although only a lower limit of the torque n/vs r 

, , • r ,r j_i where hi X ) stands for the integer part of A, and where A 

value can be inferred from the runs. . \, ' , . 

is the radial zone width. The co-orbital corotation torque 

in this situation can therefore be expressed as: 



(A.2) 



Appendix A: Finite resolution effects on the 

co-orbital corotation torque estimate 

As the radial width of the co-orbital region in the simu- 
lations presented in this work amounts to a small number 
of zone widths (typically 3 to 18), it is of interest to inves- 
tigate how bad or good is the code used here in giving an 
estimate of the co-orbital corotation torque. For this pur- 
pose, one can consider a simplified situation, such as the 
one depicted at Fig. AT. This situation neglects the grid 
curvature and describes the gas flow around the planet in 
a shearing sheet. In this situation, the corotation torque 
is given by 



pnum pnum T 

ss ss 

where: 



i=o 



— QpiA 2 ■ Q. p r p AiT, 
1 \ 3 



(A.3) 



(A.4) 



+ 



X s 

a" 



2 2 



1)^A%£ JA/+1 



is the torque due to the horseshoe fluid elements flowing 
from iA > to — iA, and where: 



— QpiA 2 ■ fl p r p AiTii 



(A.5) 
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is the torque due to the horseshoe fluid elements flowing 
from iA < to — iA. If one writes: 



<9H 

E ?; = E + iA— + 0[(iA/r p ) 2 } 



then combining Eqs. (A. 4) and (A.5) one is led to 



ay 

iM 



(A.6) 



(A.7) 



E* 3 + ( 



The relative error on the corotation torque evaluation, £ 
(T™ m - r ss )/r ss , is therefore given by: 



eon 



4 



,i=0 



(A. 



where y = x s /A is the horseshoe region half- width, ex- 
pressed in zone radial widths. The graph of the function 



£{Y) is the solid line of Fig. A.3 



The second step of this error estimate consists of eval- 
uating the discrepancy between the actual and effective 
separatrix position. Indeed, the brackets of the second line 



of Eqs. (A. 4) and (A.5) stand for the mass fraction of the 
outermost horseshoe zone included within the separatrix 
(i.e. it is also the ratio of the shaded area surface to the 



zone surface in Fig. A.l). There is no warranty however 



that the separatrix position which has to be used in this 
evaluation coincides with the one provided by the stream- 
line analysis. In other words, one has to evaluate how much 
separatrix crossing numerical effects can lead to, in order 
to infer how distant the actual and effective separatrices 
can be. One approximate way to do that is to consider 
a ID problem (along the radial dimension). Disregarding 
azimuthal advection should not lead to significant errors, 
as this latter is treated apart during the advection pro- 
cedure based on the operator splitting technique, and as 
azimuthal motion should not lead to separatrix crossing 
provided this latter is horizontal (or weakly tilted) in the 
(9, r — r p ) plane. Following the notations of Fig. |A.2| , one 
can write the total mass M_ contained below the separa- 
trix (or below any arbitrary position x s ) at time t as: 



M_ 



= A5^Sj + Ei(a,-iA) 



(A.9) 



where A is the zone spacing. This mass will be referred to 
hereafter as the inner mass. At time t + St, the inner mass 
becomes: 



M'_ = A^Ej+E^-iA) 

3=0 



(A.10) 



where the primes denote the new quantities, and where 
one has the following relationship: 



St 



(A.H) 
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Fig. A. 2. Sketch of the ID mesh used to evaluate the 
impact of numerical effects on separatrix crossing. The 
linear density E is zone-centered, while the radial velocity 
u is defined at zone interfaces (on a staggered mesh). 



where Fj is the mass flux (oriented upwards) across the 
interface between the zones J — 1 and j. One also needs 
the following relationship: 



St 



(A.12) 



which comes from the fact that the streamlines (and there- 
fore the separatrix) are found integrating the linearly in- 
terpolated velocity field. One can therefore write: 

i-l i-l 

ML = Aj2Vj+5tJ2(Fj-Fj+i) 

3=0 3=0 

St, 
-(Fi-F i+1 



A 



(x' s - iA) 



(A.13) 



The mass fluxes can be written as: 



Fj 



(A.14) 



where E* is an estimate of the linear density at the in- 
terface, the expression of which depends on the numeri- 
cal method (Stone & Norman 1992). To lowest order in 
UjSt/ A, the mass variation can be expressed, after a few 
transformations, and assuming Fq = 0, as: 



5M- 



St 



Ui(T,i - E 
«i+i(Ei - E 



)( 



i+lJ 



t + 1- 

A 



A 



(A.15) 



This mass variation necessarily corresponds to the separa- 
trix crossing mass, since Fq = 0. One can get an estimate 
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Fig. A. 3. Relative error on the co-orbital corotation 
torque in the shearing sheet, as a function of resolution 
(number of zones in the horseshoe region half- width) . The 
dotted and dashed lines are respectively for the worst case 
and a conservative case for the runs presented here (see 
text for details). 

of this mass variation assuming m — = u (a reason- 
able assumption in the case of a radial viscous drift) and 
writing Ej — S* = — ~ d r YiA/2 (in which case one 
has to assume that \u5t\ <C A, a reasonable assumption 
for the case of interest here). One is then led to: 



5M_ ~ uStd r T,(x C i - x s ) 



(A.16) 



where x\ = (i + 1/2) A is the x value of the zone center. 
One can therefore write: 



dM- 
dt 



u(x1 — x s )d r T, = x s (x^ — x s )d r Ti 



from which one can infer: 
1 



M_ = M c 



\{x\-x.ya r 'L 



(A.17) 



(A.18) 



where M°_ is the inner mass when the separatrix crosses 
the zone center. Eq. ( A.18| ) shows that the inner mass at 
the top (x s = x z c + A/2) and the bottom (x s — x l c — A/2) 
of the zone are the same, from which one can conclude 
that, although some mass can cross the separatrix as its 
sweeps the zone, there is no net effect as it sweeps the 
entire zone from one interface to the other. As a conse- 
quence, the maximum amount of mass that can cross the 
separatrix as it sweeps an arbitrary number of zones is 
M sc = d r Il ■ A 2 /8. Denoting with Sx s the error on the 
separatrix position, which is defined as: Sx s ■ £ = M SCl 
one is led to the following relation: 



A 2 d r T, 
~8 XT 



(A.19) 



The worst case is obtained for d r Y, ~ S/A, i.e. when 
the separatrix lies on the edge of a fully emptied gap, 
and when this edge is not resolved (i.e. when the density 
falls abruptly from £ to from one zone to its neighbor). 



It should be stressed that this situation is far from be- 
ing met in the simulations presented here. Even in this 
worst case, the maximal error on the separatrix position 
will be: Sx s = A/ 8, whereas in a more reasonable case 
for which d r T, will typically amount to a small fraction 
of £/£?, the maximal error on the separatrix position 
will be a fraction of (A/H) ■ A/8. The worst case is dis- 
played in Fig. A. 3 with dotted lines, while the conserva- 



tive case where Sx s = A/16, obtained for <9 r £ = E/iJ and 
H = 0.03 in the normal resolution runs, is displayed with 
dashed lines. It is worth noting that the separatrix error is 
much smaller than this conservative estimate, in particu- 
lar for small planet masses, which hardly perturb the disk 
density profile, and for which Y = x s /A is the smallest. 
One can conclude from Fig. A. 3 that the corotation torque 
is described with an accuracy better than 15 % as soon 
as x s /A > 2.3. Although this condition is not met for 
the smallest mass planets in the normal resolution runs 
presented here, the corotation torque, which on average 
for 1 < Y < 2 numerical effects tend to overestimate, is 
never found to be large enough to cancel out the differen- 
tial Lindblad torque for these objects. Therefore, one can 
conclude that even with a relatively small number of mesh 
rings involved in the co-orbital region, the code used here 
can describe with a reasonable precision the dynamics of 
the co-orbital corotation torque, and that whenever the 
planet mass is too small for the dynamics of this region to 
be correctly described (for the resolution and parameter 
set used here), the corotation torque is too small to sig- 
nificantly affect the migration anyway (one recovers the 
linear regime, see section [j] for details) . 
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